{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting Linear Models with Custom Loss Functions and Regularization in Python\n",
    "\n",
    "### by [Alex P. Miller](https://alex.miller.im/) ([@alexpmil](https://twitter.com/alexpmil))\n",
    "\n",
    "As part of a predictive model competition I participated in earlier this month, I found myself trying to accomplish a peculiar task. The challenge organizers were going to use \"mean absolute percentage error\" (MAPE) as their criterion for model evaluation. Since this is not a standard loss function built into most software, I decided to write my own code to train a model that would use the MAPE in its objective function. \n",
    "\n",
    "I started by searching through the SciKit-Learn documentation on [linear models](http://scikit-learn.org/stable/modules/classes.html#module-sklearn.linear_model) to see if the model I needed has already been developed somewhere. I thought that the `sklearn.linear_model.RidgeCV` class would accomplish what I wanted (MAPE minimization with L2 regularization), but I could not get the `scoring` argument (which supposedly lets you pass a custom loss function to the model class) to behave as I expected it to.\n",
    "\n",
    "While I highly recommend searching through existing packages to see if the model you want already exists, you should (in theory) be able to use this notebook as a template for a building linear models with an arbitrary loss function and regularization scheme."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python Code\n",
    "\n",
    "I'll be using a Jupyter Notebook (running Python 3) to build my model. If you're reading this on my website, you can find [the raw .ipynb file linked here](https://github.com/alexmill/website_notebooks/blob/master/custom-loss-function-regularization-python.ipynb); you can also run a fully-exectuable version of the notebook the the Binder platform by [clicking here](https://mybinder.org/v2/gh/alexmill/website_notebooks/master?filepath=custom-loss-function-regularization-python.ipynb).\n",
    "\n",
    "We'll start with some basic imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "from sklearn.preprocessing import StandardScaler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Your Data\n",
    "\n",
    "For the purposes of this walkthrough, I'll need to generate some raw data. Presumably, if you've found yourself here, you will want to substitute this step with one where you load your own data. \n",
    "\n",
    "I am simulating a scenario where I have 100 observations on 10 features (9 features and an intercept). The \"true\" function will simply be a linear function of these features: $y=X\\beta$. However, we want to simulate observing these data with noise. Because I'm mostly going to be focusing on the \"mean absolute percentage error\" loss function, I want my noise to be on an exponential scale, which is why I am taking exponents/logs below:\n",
    "\n",
    "$$y = e^{log(X\\beta) + \\varepsilon}, \\; \\varepsilon \\sim \\mathcal{N}(0, 0.2)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate predictors\n",
    "X_raw = np.random.random(100*9)\n",
    "X_raw = np.reshape(X_raw, (100, 9))\n",
    "\n",
    "# Standardize the predictors\n",
    "scaler = StandardScaler().fit(X_raw)\n",
    "X = scaler.transform(X_raw)\n",
    "\n",
    "# Add an intercept column to the model.\n",
    "X = np.abs(np.concatenate((np.ones((X.shape[0],1)), X), axis=1))\n",
    "\n",
    "# Define my \"true\" beta coefficients\n",
    "beta = np.array([2,6,7,3,5,7,1,2,2,8])\n",
    "\n",
    "# Y = Xb\n",
    "Y_true = np.matmul(X,beta)\n",
    "\n",
    "# Observed data with noise\n",
    "Y = Y_true*np.exp(np.random.normal(loc=0.0, scale=0.2, size=100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define your custom loss function\n",
    "\n",
    "I am mainly going to focus on the MAPE loss function in this notebook, but this is where you would substitute in your own loss function (if applicable). MAPE is defined as follows:\n",
    "\n",
    "### Mean Absolute Percentage Error (MAPE)\n",
    "\n",
    "$$ \\text{error}(\\beta) = \\frac{100}{n} \\sum_{i=1}^{n}\\left| \\frac{y_i - X_i\\beta}{y_i} \\right|$$\n",
    "\n",
    "While I won't go to into too much detail here, I ended up using a *weighted* MAPE criteria to fit the model I used in the data science competition. Given a set of sample weights $w_i$, you can define the weighted MAPE loss function using the following formula:\n",
    "\n",
    "### Weighted MAPE\n",
    "\n",
    "$$\\text{error}(\\beta) =  100 \\left( \\sum_{i=1}^N w_i \\right)^{-1} \\sum_{i=1}^N w_i \\left| \\frac{y_i - X_i\\beta}{y_i} \\right|$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Python, the MAPE can be calculated with the function below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_absolute_percentage_error(y_pred, y_true, sample_weights=None):\n",
    "    y_true = np.array(y_true)\n",
    "    y_pred = np.array(y_pred)\n",
    "    assert len(y_true) == len(y_pred)\n",
    "    \n",
    "    if np.any(y_true==0):\n",
    "        print(\"Found zeroes in y_true. MAPE undefined. Removing from set...\")\n",
    "        idx = np.where(y_true==0)\n",
    "        y_true = np.delete(y_true, idx)\n",
    "        y_pred = np.delete(y_pred, idx)\n",
    "        if type(sample_weights) != type(None):\n",
    "            sample_weights = np.array(sample_weights)\n",
    "            sample_weights = np.delete(sample_weights, idx)\n",
    "        \n",
    "    if type(sample_weights) == type(None):\n",
    "        return(np.mean(np.abs((y_true - y_pred) / y_true)) * 100)\n",
    "    else:\n",
    "        sample_weights = np.array(sample_weights)\n",
    "        assert len(sample_weights) == len(y_true)\n",
    "        return(100/sum(sample_weights)*np.dot(\n",
    "                sample_weights, (np.abs((y_true - y_pred) / y_true))\n",
    "        ))\n",
    "    \n",
    "loss_function = mean_absolute_percentage_error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting a simple linear model with custom loss function\n",
    "\n",
    "You may know that the traditional method for fitting linear models, ordinary least squares, has a nice analytic solution. This means that the \"optimal\" model parameters that minimize the squared error of the model, can be calculated directly from the input data:\n",
    "\n",
    "$$ \\hat\\beta = \\arg\\min_\\beta \\frac{1}{n} \\sum_{i=1}^n (y_i - X_i\\beta)^2 =  (X^\\mathrm{T}X)^{-1}X^\\mathrm{T}y $$\n",
    "\n",
    "However, with an arbitrary loss function, there is no guarantee that finding the optimal parameters can be done so easily. To keep this notebook as generalizable as possible, I'm going to be minimizing our custom loss functions using numerical optimization techniques (similar to the \"solver\" functionality in Excel). In general terms, the $\\beta$ we want to fit can be found as the solution to the following equation (where I've subsituted in the MAPE for the error function in the last line):\n",
    "\n",
    "$$ \\hat\\beta = \\arg\\min_\\beta \\; \\text{error}(\\beta) = \\arg\\min_\\beta \\frac{100}{n} \\sum_{i=1}^{n}\\left| \\frac{y_i - X_i\\beta}{y_i} \\right| $$\n",
    "\n",
    "Essentially we want to search over the space of all $\\beta$ values and find the value that minimizes our chosen error function. To get a flavor for what this looks like in Python, I'll fit a simple MAPE model below, using the `minimize` function from SciPy. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 9.08252394  5.54995839  8.75233095  1.1883712   3.29482497  5.03886496\n",
      " -0.22556182  0.38830739  3.15524308  5.24599191]\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "def objective_function(beta, X, Y):\n",
    "    error = loss_function(np.matmul(X,beta), Y)\n",
    "    return(error)\n",
    "\n",
    "# You must provide a starting point at which to initialize\n",
    "# the parameter search space\n",
    "beta_init = np.array([1]*X.shape[1])\n",
    "result = minimize(objective_function, beta_init, args=(X,Y),\n",
    "                  method='BFGS', options={'maxiter': 500})\n",
    "\n",
    "# The optimal values for the input parameters are stored\n",
    "# in result.x\n",
    "beta_hat = result.x\n",
    "print(beta_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare the esimated betas to the true model betas that we initialized at the beginning of this notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>true_beta</th>\n",
       "      <th>estimated_beta</th>\n",
       "      <th>error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2</td>\n",
       "      <td>9.082524</td>\n",
       "      <td>-7.082524</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>6</td>\n",
       "      <td>5.549958</td>\n",
       "      <td>0.450042</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>7</td>\n",
       "      <td>8.752331</td>\n",
       "      <td>-1.752331</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>1.188371</td>\n",
       "      <td>1.811629</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>3.294825</td>\n",
       "      <td>1.705175</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>7</td>\n",
       "      <td>5.038865</td>\n",
       "      <td>1.961135</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1</td>\n",
       "      <td>-0.225562</td>\n",
       "      <td>1.225562</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>2</td>\n",
       "      <td>0.388307</td>\n",
       "      <td>1.611693</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2</td>\n",
       "      <td>3.155243</td>\n",
       "      <td>-1.155243</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>8</td>\n",
       "      <td>5.245992</td>\n",
       "      <td>2.754008</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   true_beta  estimated_beta     error\n",
       "0          2        9.082524 -7.082524\n",
       "1          6        5.549958  0.450042\n",
       "2          7        8.752331 -1.752331\n",
       "3          3        1.188371  1.811629\n",
       "4          5        3.294825  1.705175\n",
       "5          7        5.038865  1.961135\n",
       "6          1       -0.225562  1.225562\n",
       "7          2        0.388307  1.611693\n",
       "8          2        3.155243 -1.155243\n",
       "9          8        5.245992  2.754008"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame({\n",
    "    \"true_beta\": beta, \n",
    "    \"estimated_beta\": beta_hat,\n",
    "    \"error\": beta-beta_hat\n",
    "})[[\"true_beta\", \"estimated_beta\", \"error\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's obviously not perfect, but we can see that our estimated values are at least in the ballpark from the true values. We can also calculate the final MAPE of our estimated model using our loss function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14.354033248368872"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loss_function(np.matmul(X,beta_hat), Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Incorporating Regularization into Model Fitting\n",
    "\n",
    "The process described above fits a simple linear model to the data provided by directly minimizing the a custom loss function (MAPE, in this case). However, in many machine learning problems, you will want to [regularize](https://en.wikipedia.org/wiki/Regularization_(mathematics)) your model parameters to prevent overfitting. In this notebook, I'm going to walk through the process of incorporating L2 regularization, which amounts to penalizing your model's parameters by the square of their magnitude. \n",
    "\n",
    "In precise terms, rather than minimizing our loss function directly, we will augment our loss function by adding a squared penalty term on our model's coefficients. With L2 regularization, our new loss function becomes:\n",
    "\n",
    "$$ L(\\beta) = \\frac{100}{N} \\sum_{i=1}^N \\left| \\frac{y_i - X_i\\beta}{y_i} \\right| + \\lambda \\sum_{k=1}^K \\beta_k^2 $$\n",
    "\n",
    "Or, in the case that sample weights are provided:\n",
    "\n",
    "$$ L(\\beta) = 100 \\left( \\sum_{i=1}^N w_i \\right)^{-1} \\sum_{i=1}^N w_i \\left| \\frac{y_i - X_i\\beta}{y_i} \\right| + \\lambda \\sum_{k=1}^K \\beta_k^2 $$\n",
    "\n",
    "For now, we will assume that the $\\lambda$ coefficient (the regularization parameter) is already known. However, later we will use cross validation to find the optimal $\\lambda$ value for our data.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since our  model is getting a little more complicated, I'm going to define a Python class with a very similar attribute and method scheme as those found in SciKit-Learn (e.g., `sklearn.linear_model.Lasso` or `sklearn.ensemble.RandomForestRegressor`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CustomLinearModel:\n",
    "    \"\"\"\n",
    "    Linear model: Y = XB, fit by minimizing the provided loss_function\n",
    "    with L2 regularization\n",
    "    \"\"\"\n",
    "    def __init__(self, loss_function=mean_absolute_percentage_error, \n",
    "                 X=None, Y=None, sample_weights=None, beta_init=None, \n",
    "                 regularization=0.00012):\n",
    "        self.regularization = regularization\n",
    "        self.beta = None\n",
    "        self.loss_function = loss_function\n",
    "        self.sample_weights = sample_weights\n",
    "        self.beta_init = beta_init\n",
    "        \n",
    "        self.X = X\n",
    "        self.Y = Y\n",
    "            \n",
    "    \n",
    "    def predict(self, X):\n",
    "        prediction = np.matmul(X, self.beta)\n",
    "        return(prediction)\n",
    "\n",
    "    def model_error(self):\n",
    "        error = self.loss_function(\n",
    "            self.predict(self.X), self.Y, sample_weights=self.sample_weights\n",
    "        )\n",
    "        return(error)\n",
    "    \n",
    "    def l2_regularized_loss(self, beta):\n",
    "        self.beta = beta\n",
    "        return(self.model_error() + \\\n",
    "               sum(self.regularization*np.array(self.beta)**2))\n",
    "    \n",
    "    def fit(self, maxiter=250):        \n",
    "        # Initialize beta estimates (you may need to normalize\n",
    "        # your data and choose smarter initialization values\n",
    "        # depending on the shape of your loss function)\n",
    "        if type(self.beta_init)==type(None):\n",
    "            # set beta_init = 1 for every feature\n",
    "            self.beta_init = np.array([1]*self.X.shape[1])\n",
    "        else: \n",
    "            # Use provided initial values\n",
    "            pass\n",
    "            \n",
    "        if self.beta!=None and all(self.beta_init == self.beta):\n",
    "            print(\"Model already fit once; continuing fit with more itrations.\")\n",
    "            \n",
    "        res = minimize(self.l2_regularized_loss, self.beta_init,\n",
    "                       method='BFGS', options={'maxiter': 500})\n",
    "        self.beta = res.x\n",
    "        self.beta_init = self.beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 8.70454035,  5.56955027,  8.82671937,  1.10660836,  3.36271348,\n",
       "        5.16710648, -0.08675964,  0.4776243 ,  3.12646051,  5.28643399])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l2_mape_model = CustomLinearModel(\n",
    "    loss_function=mean_absolute_percentage_error,\n",
    "    X=X, Y=Y, regularization=0.00012\n",
    ")\n",
    "l2_mape_model.fit()\n",
    "l2_mape_model.beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just to confirm that our regularization did work, let's make sure that the estimated betas found with regularization are different from those found without regularization (which we calculated earlier):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>true_beta</th>\n",
       "      <th>estimated_beta</th>\n",
       "      <th>regularized_beta</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2</td>\n",
       "      <td>9.082524</td>\n",
       "      <td>8.704540</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>6</td>\n",
       "      <td>5.549958</td>\n",
       "      <td>5.569550</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>7</td>\n",
       "      <td>8.752331</td>\n",
       "      <td>8.826719</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>1.188371</td>\n",
       "      <td>1.106608</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>3.294825</td>\n",
       "      <td>3.362713</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>7</td>\n",
       "      <td>5.038865</td>\n",
       "      <td>5.167106</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1</td>\n",
       "      <td>-0.225562</td>\n",
       "      <td>-0.086760</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>2</td>\n",
       "      <td>0.388307</td>\n",
       "      <td>0.477624</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2</td>\n",
       "      <td>3.155243</td>\n",
       "      <td>3.126461</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>8</td>\n",
       "      <td>5.245992</td>\n",
       "      <td>5.286434</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   true_beta  estimated_beta  regularized_beta\n",
       "0          2        9.082524          8.704540\n",
       "1          6        5.549958          5.569550\n",
       "2          7        8.752331          8.826719\n",
       "3          3        1.188371          1.106608\n",
       "4          5        3.294825          3.362713\n",
       "5          7        5.038865          5.167106\n",
       "6          1       -0.225562         -0.086760\n",
       "7          2        0.388307          0.477624\n",
       "8          2        3.155243          3.126461\n",
       "9          8        5.245992          5.286434"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame({\n",
    "    \"true_beta\": beta, \n",
    "    \"estimated_beta\": beta_hat,\n",
    "    \"regularized_beta\": l2_mape_model.beta\n",
    "})[[\"true_beta\", \"estimated_beta\", \"regularized_beta\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since our regularization parameter is so small, we can see that it didn't affect our coefficient estimates dramatically. But the fact that the betas are different between the two models indicates that our regularization does seem to be working."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just to make sure things are in the realm of common sense, it's never a bad idea to plot your predicted Y against our observed Y."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7e09da60a8d0>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGy1JREFUeJzt3X+MXPdZ7/H348202aQ069w6lrNpcG8VuVyIsJtV6ZUlRBPalLZK9qZKf0hF5ipS+IN7VXIr0w1CkKKiLAQo/IUwBWTdBm5Mk2zSFhoiOwU1ormsuw65uY4VAWlg7WtvIQsN3tKN/dw/9owzO3vOzJmZ8+N7vvN5SdbuzM7uPHN2/cz3PN/n+z3m7oiISPNtqzsAEREphhK6iEgklNBFRCKhhC4iEgkldBGRSCihi4hEQgldRCQSSugiIpFQQhcRicRlVT7ZW97yFt+9e3eVTyki0njHjx//trvv6Pe4ShP67t27WVxcrPIpRUQaz8y+ledxKrmIiERCCV1EJBJK6CIikVBCFxGJhBK6iEgkKu1yEZHxsrC0zANPnOL06hrXTk1y8NY9zO6brjusaCmhi0gpFpaWufeR51hbvwDA8uoa9z7yHICSeklUchGRUjzwxKlLybxtbf0CDzxxqqaI4qeELiKlOL26NtD9MjoldBEpxbVTkwPdL6NTQheRUhy8dQ+TrYlN9022Jjh4656aIoqfJkVFpBTtiU91uVRHCV1ESjO7b1oJvEIquYiIREIJXUQkEkroIiKRUEIXEYmEErqISCSU0EVEIqGELiISCSV0EZFIKKGLiESib0I3sz1mdqLj37+a2c+a2dVm9qSZvZh83F5FwCIikq5vQnf3U+6+1933AjcB54FHgTngqLvfABxNbouISE0GLbncAvytu38LuB04nNx/GJgtMjARERnMoAn9Y8AfJ5/vdPczAMnHa4oMTEREBpM7oZvZG4DbgD8Z5AnM7G4zWzSzxZWVlUHjExGRnAYZof8E8E13P5vcPmtmuwCSj+fSvsndD7n7jLvP7NixY7RoRUQk0yAJ/eO8Xm4BeBw4kHx+AHisqKBERGRwuRK6mV0BvBd4pOPueeC9ZvZi8rX54sMTEZG8cl2xyN3PA/+h675/YqPrRUREAqCVoiIikVBCFxGJhBK6iEgklNBFRCKhhC4iEgkldBGRSCihi4hEQgldRCQSSugiIpFQQhcRiUSupf8iIk2ysLTMA0+c4vTqGtdOTXLw1j3M7puuO6zSKaGLSFQWlpa595HnWFu/AMDy6hr3PvIcQPRJXSUXEYnKA0+cupTM29bWL/DAE6dqiqg6SugiEpXTq2sD3R8TJXQRicq1U5MD3R8TJXSRCCwsLbN//hhvm/sK++ePsbC0XHdItTl46x4mWxOb7ptsTXDw1j01RVQdTYqKNNw4TwKmab9mdbmISOP0mgQchySWZnbf9Fi+dpVcRBpunCcBZTMldJGGG+dJQNlMCV2k4cZ5ElA2Uw1dpOHGeRJQNlNCF4nAuE4CymYquYiIRCJXQjezKTP7opm9YGYnzew/m9nVZvakmb2YfNxedrAiIpIt7wj9t4Gvuvs7gB8GTgJzwFF3vwE4mtwWEZGa9K2hm9mbgR8FfgrA3b8HfM/Mbgd+LHnYYeBrwKfLCFJEwjeue5CHJM8I/T8CK8AfmtmSmX3ezK4Edrr7GYDk4zUlxikiAWtvP7C8uobz+vYD47ynTB3yJPTLgHcCv+Pu+4B/Y4DyipndbWaLZra4srIyZJgiErJx3oM8JHkS+j8C/+juzyS3v8hGgj9rZrsAko/n0r7Z3Q+5+4y7z+zYsaOImEUkMNp+IAx9a+ju/v/M7B/MbI+7nwJuAf5v8u8AMJ98fKzUSEUCFXrtuIr4rp2aZDkleWv7gWrlXVj034EHzewNwN8B/5WN0f0RM7sLeBm4s5wQRcIV+ta1VcV38NY9m54HtP1AHXIldHc/AcykfOmWYsMRaZbQt66tKj5tP5Cu6rM3Lf0XGUHoteMq49P2A5vVcfampf8iIwh969rQ44tZHZ0/SugiIwh969rQ44tZHWdvKrmIjCD02nHo8cWsjs4fc/fSfni3mZkZX1xcrOz5RETq0l1Dh42zo/vvuHHgN1QzO+7uaY0pm2iELiJSgjrOjpTQRaQUoS+4qkLVnT9K6CIlG8fEFvqCq1ipy0WkROO6C6E266qHRugiJQppJWmVZwqhL7iKlRK6SIlCSWxVl0CatllXLGUxlVxEShTKSs2qSyBNWtAUU1lMCV2kRKEktqrPFGb3TXP/HTcyPTWJAdNTk0P1X1chpnq/Si4iJQplpWaRJZC85YmmbNYVSlmsCEroIiULIbEVtV95jO2ITav396KSi8gYKKoEElN5oi2UslgRNEIXGRNFnCnEVJ5oC6UsVgQldBHJLabyRKcQymJFUMlFRHKLqTwRI43QRSS3mMoTMVJCFxlRLKsM84qlPBEjJXSJWtnJNsY2PmkuJXSJVhXJNqTNt6R+dZ+tKaFLtKpItjG28clwiTmEs7VcCd3MXgK+A1wAXnP3GTO7GngI2A28BHzE3V8pJ0yRwVWRbGNt40tT9+izKsMm5hDO1gZpW3yPu+/tuFDpHHDU3W8Ajia3RYJRxU6H49LGF9OOhP0Muxo2hLO1UfrQbwcOJ58fBmZHD0ekv4WlZfbPH+Ntc19h//yxzKRSRbKtalfBvK+5LDEu+c8ybGIOYavkvDV0B/7czBz4XXc/BOx09zMA7n7GzK5J+0Yzuxu4G+D6668vIGQZZ4OcDlfVM112G18ItdkQRp9VGbaMVtQGaKPIm9D3u/vpJGk/aWYv5H2CJPkfApiZmfEhYhS5ZNA6ZQw90yHUZkOcKyirpj9sYg5h0VWuhO7up5OP58zsUeBdwFkz25WMzncB50qMUwQYr5FiWwivOYTRZ6eiz1q63xw+fNM0T72wMnBirnsA0Tehm9mVwDZ3/07y+fuAXwYeBw4A88nHx8oMVATCHCmWLYTXHMLos1ORZy1pbw4PH18O9gpLveQZoe8EHjWz9uP/yN2/amZ/DRwxs7uAl4E7ywtTZENoI8UqhPKa6x59diryrCWEklZR+iZ0d/874IdT7v8n4JYyghLJEtpIsQrj+Jr7KfKsJYSSVlG0UlQaJ6SRYlXG8TX3UuRZSwglraJoP3QRaZwi+/9jWhymEbpIQ4zL0vu8ijpriamkpYQu0gAhLC6KWSwlLZVcRBpgnJbey/CU0EUaIKZODCmPSi4iDRB6J4bq+2HQCF2kAULuxGjS1rp171pZNiV0kQaoapveYTSlvt+kN55hqeQi0hChdmI0pb4f0xL/LBqhi8hIQriwQx5NeeMZhRK6iIwk5Pp+p6a88YxCCV1ERhJyfb9TU954RqEausiI1LI3eH2/jmMW0xL/LEroIiPQkvzB1XnMQp1YLooSusgIQuic6DfaDe0MIoRjFisldJER1N050W+0G+IZRN3HLGaaFJXSxbw6r+7OiX6LekJc9FP3MYuZErqUqkmr84Z54ymrcyJvLP1GuyGOhseh26QuKrmMgTprqEXVS8t+DcOWJorqnOh8fVNXtHj1u6+xftH7xtJv064QN/Uah26TuiihR67uGmoRI8QqXsMobzyjdk50v75Xzq9vecza+gU+86XntzxPv2trFnntzSLF3m1SF5VcIld3DbWIemkVr6HO0kTa60vzyvl1dneVYPot6pndN82Hb5pmwgyACTM+fJOSaaw0Qo9c3TXUIkaIVbyGOksTg76O7jOUXqPdhaVlHj6+zAXfKN9ccOfh48vMfP/VQZW8pBi5R+hmNmFmS2b25eT21Wb2pJm9mHzcXl6YMqy6OwqKWBZe1GvoNdFY50TdML+LvGcoRZzdNGlie9wNUnL5JHCy4/YccNTdbwCOJrclMCF0FMzum+bpuZv5+/kP8vTczQOP7Ip4Df2SUvcbz/YrWrzxsm3c89CJ0lst015fHnlG9kWc3dRdtpP8ciV0M7sO+CDw+Y67bwcOJ58fBmaLDU1G0R6N3vPQCS5vbWNqshX0xkm9FDHKz5OU2m88n/voXr67fpHVtfVKRqTt1zc12Rro+/KM7Is4u6m7bCf55a2h/xbwc8D3ddy3093PALj7GTO7pujgZDhpXROTrQk+99G9jUrknUbtihgkKdWxNH123zQPPHGK1bWtHS4ABnjH7bxnKEXMYVw12UqNSwuBwtN3hG5mHwLOufvxYZ7AzO42s0UzW1xZWRnmR8iAdIq81SAj1bpGpFk/34DPfXTvUGcoaWc377z+Kj515Fl2z32Ft9/7p/zCwnOZ37+wtMy/fe+1Lfe3tlntrY+yVZ4R+n7gNjP7AHA58GYz+wJw1sx2JaPzXcC5tG9290PAIYCZmRlPe4wUS6fIWw0yUq2r46XX845yhtL5vb+w8Bxf+MbLl752wf3S7c/O3rjlex944hTrF7b+t33T5Zc19mwvZn1H6O5+r7tf5+67gY8Bx9z9E8DjwIHkYQeAx0qLUgZSd2dLEYre/2WQOnxdE8lVPO8fP/MPA92fNQh45fx6dPvyxGCUPvR54IiZ3QW8DNxZTEgyqlBXB+Y1yMrQPP3R3Y/pN5dQ19L0Kp633Y+e9/6sswYIY+dG2cw84xdZhpmZGV9cXKzs+cZZkxeC7J8/lppEpqcmeXru5ku3uxM/bLxxdY688zxmUE0+tm+/909Tk/eEGX97/we23J92/Lp1/17yaPIxrIOZHXf3mX6P00rRSDV5r4y8cwBZk7+fOvIs8HrnSJEdK3XvjTOqj//IWzfV0DvvT9N51pA1Uh90bqbpxzBk2stFgpN3DiArkVxwv9Q3XvQEcdM7iD47eyOfePf1m/Z2+cS7r0+dEG1r9+dPFzQ30/RjGDKN0CU4eecAetV32wlimI6VXuWAGDqIPjt7Y88EnqWouZkYjmGoNEKX4OTtSOm3ZP706trAnSP9tgiIoYNoWEWs2IXxPoZl0whdgpRnDqD99U8deTZ1oq/dvw3pnSNpI/F+NfemdhAtLC3zmS89f2mv9anJFvfd9oMDJ+Mi5maaegybQAldGq2dXHoliLQklDUxl9XN0S4HNPFqOwtLyxz84rObFgitrq1z8E9enzyuUhOPYVMooUvjdXdiTJhtmmRLSxRZI/EJs8zRfufzNSn5ZK32XL/ope5P00vTjmFTqIYuUWiXQyZbE5cS8vLqGvc8dCJ1r5JeHTJ1bzdctF6TjZqIjIsSukQjbdTtwIPfeHnLEvWsCbj2RN+oE38h6TXZqInIuKjkIj01aUVf1mjTYUtpodfEXNXlgLKP8cFb92ypoYN2TIyRErpkatqKvl596d3JPpSJuSqOcfvnFNHlImHTXi6SKe+eKqFYWFrmnodOkPYXHWrMTTvGUg/t5SIjC2FF3yDliNl90yx+65958BsvD3V1nzqEcIwlHpoUlUx1r+gb5mrzn529ceir+9Sh7mMscdEIXTLVvaJv2J0Sm9TjXPcxlrhohC6Zitq7Y1jjUI7oPsZTky0ub23jnodO6IpAMjCN0BuoylbCOke7dV3bM48ifwftY9y0riIJjxJ6w4T4n75746crWtt4w2UT/Mvaeu7LwqU9pqpyxKDJudfvAIZvhSz6YhwyfpTQGya0//RpGz+dX7/I+fWLQPobTt43pSp6xQd9g1xYWk7d3XFt/QL3Pf48//7axaHfbMehxCTlUkKvSPco8D3v2MFTL6wMnKhC+0+ftfFTp1EuC1d2yWeQWNrJP+uCyqtr61vuG+TNNuQSkzSDEnoF0kaBndd1zBrJpZUCQvtPn7Uys1v7snAQ1pvSILGkJf9RnqObOl5kVOpyqUCeRNB9TcWsHuz3vGNHMLsBLiwtYwM8vj1SzxrP1/GmNEgfeK/EPNmaYPsVrYGeo1vdXUXSfBqhVyDvCK3zcVmlgKdeWOH+O27MfQWeopJB1tV9Bt04IqtcUdeb0iCj4qyzowkz7r9j4xqdo46wm9RDL+FRQq9Ar02juh/X1qsUMMgVeNoGTfSdCXzqihavfvc11i++vs94r6v7ABjkTvbTFW+M1f3m9OGbpnPNZ2Ql/+5RdN0bfsn46pvQzexy4C+BNyaP/6K7/5KZXQ08BOwGXgI+4u6vlBdqc6Ulgm7dI7mrJlupk2xXTaaf1meN6D/zpef57vpgnRfdbw7tdsTun511dZ92gu73mmEj8Ve5CVXaG9/Dx5dzlTbydN1ohC11yjNC/3fgZnd/1cxawNfN7M+AO4Cj7j5vZnPAHPDpEmOtXFEljLRE0K/LxTKK01n3Z43os5Jxr86LvJN/7av7ZO0p3v5Zp1fX2Jbj0m5VGLXtszNht/8+7nnohEbjEoS+Cd039td9NbnZSv45cDvwY8n9h4GvEVFCL3oBz6Ajt9WURNzr/rxlnbbl1TX2zx9LfUPJW/Of7qilp/2c7uQXQgdHUR02IS7wEsnV5WJmE2Z2AjgHPOnuzwA73f0MQPLxmvLCrF6vkVwVBt2Fr309zU6TrQmmMko0Bpm7GOYZNXeOxJ+eu5m/n/8gT8/d3HNr2xA6OIra3bDuvw+RNLkSurtfcPe9wHXAu8zsh/I+gZndbWaLZra4srIybJyVq7tXOitBZ41osxLmfbf94JafA1snLDuTUdpztyaMqcnWSMk4b/Iv06DHNUvdfx8iaQbqcnH3VTP7GvB+4KyZ7XL3M2a2i43Re9r3HAIOwcYVi0aMtzJ1L+AZZtl7Vlkn7aIPadrJKJTLs5WhqNdW99+HSJo8XS47gPUkmU8CPw78KvA4cACYTz4+VmagVQth1V5RHRNPvbCSq4WwMxnF3K1RxGsL4e9DpFueEfou4LCZTbBRojni7l82s78CjpjZXcDLwJ0lxlm5mEapecoAMSajMhdaxfT3IfHQRaLHQNaFiCfMuOgeZTLK6qrRUnppIl0kWi7Ju8IxJqFtMyxSBSX0MTCO5QF1ocg4UkIfEzFPcqZRF4qMIyV0iUbnJOhVky1aE7bp4hsxTvyKdIoioVd50WQZXRm/r+5J0NW1dVrbjO1XtFg9n31tU5GYND6ha0+NZinr95U2Cbp+0bniDZex9IvvGz5gkQZp/BWLtKfGVgtLy+yfP8bb5r7C/vljl/ZoCUFZvy9NgopEkND1H3mzrEvXhZLUy/p9FbXplkiTNT6h6z/yZqGfsZTx+1pYWub8917bcr8mQWXcND6hF7V7XizqOGMZpMRT9O+rfUbSfSGPqclW1AunRNI0flJ0HBfN9FJ1//Wgk5xF/76yrq505RsvG9u/ARlfjU/oMH6LZnqpehfAYZbYF/n70hyKyOsaX3KRzaq+MlDdCVVzKCKvi2KELptVecZS9xJ77Usu8jqN0GUkdU9Kh3KtUpEQaIQeqaq2QwhhUlpzKCIblNAjVPV2CEqoImFQySVCoS8uEpFyKKFHqO7OExGphxJ6hNTKJzKelNAjVHfniYjUQ5OiEQqh86RqusiJiBJ6tMap80QXORHZoIQeGI00BzfMfjIiMepbQzezt5rZU2Z20syeN7NPJvdfbWZPmtmLycft5Ycbt9AvThEqdfWIbMgzKfoa8Cl3/wHg3cDPmNl/AuaAo+5+A3A0uS0jUP/4cNTVI7Khb0J39zPu/s3k8+8AJ4Fp4HbgcPKww8BsWUGOC400h6OuHpENA7UtmtluYB/wDLDT3c/ARtIHrsn4nrvNbNHMFldWVkaLNnIaaQ5HG3SJbDB3z/dAszcBfwH8irs/Ymar7j7V8fVX3L1nHX1mZsYXFxdHCjhm3d0asDHSVHISGW9mdtzdZ/o9LleXi5m1gIeBB939keTus2a2y93PmNku4Nzw4Q4nto6QcewfF5Hi9E3oZmbA7wMn3f03O770OHAAmE8+PlZKhBli7T0ep/5xESlWnhr6fuAngZvN7ETy7wNsJPL3mtmLwHuT25VRR4iIyGZ9R+ju/nXAMr58S7Hh5BdLR0hsZSMRqU9jN+eKoSNEC4lEpEiNTegx9B6rbCQiRWrsXi4xdITEUjYSkTA0NqFD8ztCrp2aZDkleTepbCQi4WhsySUGMZSNRCQcjR6hN10MZSMRCYcSes2aXjYSkXCo5CIiEongR+haeCMikk/QCT3W/VpERMoQdMlFC29ERPILOqFr4Y2ISH5BJ/QY9msREalK0AldC29ERPILelJUC29ERPILOqGDFt6IiOQVdMlFRETyU0IXEYmEErqISCSU0EVEIqGELiISCXP36p7MbAX4Vp+HvQX4dgXhlEXx10vx16vJ8Ycc+/e7+45+D6o0oedhZovuPlN3HMNS/PVS/PVqcvxNjr1NJRcRkUgooYuIRCLEhH6o7gBGpPjrpfjr1eT4mxw7EGANXUREhhPiCF1ERIZQa0I3s7ea2VNmdtLMnjezTyb3X21mT5rZi8nH7XXGmaZH7PeZ2bKZnUj+faDuWNOY2eVm9r/N7Nkk/s8k9wd/7KFn/I04/m1mNmFmS2b25eR2I45/W0r8jTn+ZvaSmT2XxLmY3Neo49+t1pKLme0Cdrn7N83s+4DjwCzwU8A/u/u8mc0B293907UFmqJH7B8BXnX3X681wD7MzIAr3f1VM2sBXwc+CdxB4Mceesb/fhpw/NvM7H8AM8Cb3f1DZvZrNOD4t6XEfx8NOf5m9hIw4+7f7rivUce/W60jdHc/4+7fTD7/DnASmAZuBw4nDzvMRqIMSo/YG8E3vJrcbCX/nAYce+gZf2OY2XXAB4HPd9zdiOMPmfE3XWOOf5pgauhmthvYBzwD7HT3M7CROIFr6ousv67YAf6bmf2Nmf1ByKdsyenyCeAc8KS7N+rYZ8QPDTn+wG8BPwdc7LivMcef9PihOcffgT83s+NmdndyX5OO/xZBJHQzexPwMPCz7v6vdccziJTYfwd4O7AXOAP8Ro3h9eTuF9x9L3Ad8C4z+6G6YxpERvyNOP5m9iHgnLsfrzuWYfSIvxHHP7Hf3d8J/ATwM2b2o3UHNKraE3pS/3wYeNDdH0nuPpvUqNu16nN1xddLWuzufjZJNBeB3wPeVWeMebj7KvA1NurPjTj2nTrjb9Dx3w/cltRx/xdws5l9geYc/9T4G3T8cffTycdzwKNsxNqU45+q7i4XA34fOOnuv9nxpceBA8nnB4DHqo6tn6zY238Mif8C/J+qY8vDzHaY2VTy+STw48ALNODYQ3b8TTn+7n6vu1/n7ruBjwHH3P0TNOT4Z8XflONvZlcmzQyY2ZXA+9iItRHHP0vd1xTdD/wk8FxSCwX4eWAeOGJmdwEvA3fWFF8vWbF/3Mz2slGfewn46XrC62sXcNjMJth4Yz/i7l82s78i/GMP2fH/z4Yc/yxN+Nvv5dcacvx3Ao9ujMu4DPgjd/+qmf01DT7+WikqIhKJ2mvoIiJSDCV0EZFIKKGLiERCCV1EJBJK6CIikVBCFxGJhBK6iEgklNBFRCLx/wH++cyybkIW8AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Predicted Y vs. observed Y\n",
    "plt.scatter(l2_mape_model.predict(X), Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Important Caveat: Standardize Your Predictors\n",
    "> In most applications, your features will be measured on many different scales; however you'll notice in the loss function described above, each $\\beta_k$ parameter is being penalized by the same amount ($\\lambda$). Best practice when using L2 regularization is to **standardize your feature matrix** (subtract the mean off of each column and divide the result by the column standard deviation). This will ensure that all features are on approximately the same scale and that the regularization parameter has an equal impact on all $\\beta_k$ coefficients.\n",
    "\n",
    "> I standardized my data at the very beginning of this notebook, but typically you will need to work standardization into your data pipeline. Use `sklearn.preprocessing.StandardScaler` and keep track of your intercept when going through this process!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross Validation to Identify Optimal Regularization Parameter\n",
    "\n",
    "At this point, we have a model class that will find the optimal beta coefficients to minimize the loss function described above with a given regularization parameter. Of course, your regularization parameter $\\lambda$ will not typically fall from the sky. Below I've included some code that uses cross validation to find the optimal $\\lambda$, among the set of candidates provided by the user."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import KFold\n",
    "\n",
    "# Used to cross-validate models and identify optimal lambda\n",
    "class CustomCrossValidator:\n",
    "    \n",
    "    \"\"\"\n",
    "    Cross validates arbitrary model using MAPE criterion on\n",
    "    list of lambdas.\n",
    "    \"\"\"\n",
    "    def __init__(self, X, Y, ModelClass,\n",
    "                 sample_weights=None,\n",
    "                 loss_function=mean_absolute_percentage_error):\n",
    "        \n",
    "        self.X = X\n",
    "        self.Y = Y\n",
    "        self.ModelClass = ModelClass\n",
    "        self.loss_function = loss_function\n",
    "        self.sample_weights = sample_weights\n",
    "    \n",
    "    def cross_validate(self, lambdas, num_folds=10):\n",
    "        \"\"\"\n",
    "        lambdas: set of regularization parameters to try\n",
    "        num_folds: number of folds to cross-validate against\n",
    "        \"\"\"\n",
    "        \n",
    "        self.lambdas = lambdas\n",
    "        self.cv_scores = []\n",
    "        X = self.X\n",
    "        Y = self.Y \n",
    "        \n",
    "        # Beta values are not likely to differ dramatically\n",
    "        # between differnt folds. Keeping track of the estimated\n",
    "        # beta coefficients and passing them as starting values\n",
    "        # to the .fit() operator on our model class can significantly\n",
    "        # lower the time it takes for the minimize() function to run\n",
    "        beta_init = None\n",
    "        \n",
    "        for lam in self.lambdas:\n",
    "            print(\"Lambda: {}\".format(lam))\n",
    "            \n",
    "            # Split data into training/holdout sets\n",
    "            kf = KFold(n_splits=num_folds, shuffle=True)\n",
    "            kf.get_n_splits(X)\n",
    "            \n",
    "            # Keep track of the error for each holdout fold\n",
    "            k_fold_scores = []\n",
    "            \n",
    "            # Iterate over folds, using k-1 folds for training\n",
    "            # and the k-th fold for validation\n",
    "            f = 1\n",
    "            for train_index, test_index in kf.split(X):\n",
    "                # Training data\n",
    "                CV_X = X[train_index,:]\n",
    "                CV_Y = Y[train_index]\n",
    "                CV_weights = None\n",
    "                if type(self.sample_weights) != type(None):\n",
    "                    CV_weights = self.sample_weights[train_index]\n",
    "                \n",
    "                # Holdout data\n",
    "                holdout_X = X[test_index,:]\n",
    "                holdout_Y = Y[test_index]\n",
    "                holdout_weights = None\n",
    "                if type(self.sample_weights) != type(None):\n",
    "                    holdout_weights = self.sample_weights[test_index]\n",
    "                \n",
    "                # Fit model to training sample\n",
    "                lambda_fold_model = self.ModelClass(\n",
    "                    regularization=lam,\n",
    "                    X=CV_X,\n",
    "                    Y=CV_Y,\n",
    "                    sample_weights=CV_weights,\n",
    "                    beta_init=beta_init,\n",
    "                    loss_function=self.loss_function\n",
    "                )\n",
    "                lambda_fold_model.fit()\n",
    "                \n",
    "                # Extract beta values to pass as beta_init \n",
    "                # to speed up estimation of the next fold\n",
    "                beta_init = lambda_fold_model.beta\n",
    "                \n",
    "                # Calculate holdout error\n",
    "                fold_preds = lambda_fold_model.predict(holdout_X)\n",
    "                fold_mape = mean_absolute_percentage_error(\n",
    "                    holdout_Y, fold_preds, sample_weights=holdout_weights\n",
    "                )\n",
    "                k_fold_scores.append(fold_mape)\n",
    "                print(\"Fold: {}. Error: {}\".format( f, fold_mape))\n",
    "                f += 1\n",
    "            \n",
    "            # Error associated with each lambda is the average\n",
    "            # of the errors across the k folds\n",
    "            lambda_scores = np.mean(k_fold_scores)\n",
    "            print(\"LAMBDA AVERAGE: {}\".format(lambda_scores))\n",
    "            self.cv_scores.append(lambda_scores)\n",
    "        \n",
    "        # Optimal lambda is that which minimizes the cross-validation error\n",
    "        self.lambda_star_index = np.argmin(self.cv_scores)\n",
    "        self.lambda_star = self.lambdas[self.lambda_star_index]\n",
    "        print(\"\\n\\n**OPTIMAL LAMBDA: {}**\".format(self.lambda_star))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lambda: 1\n",
      "Fold: 1. Error: 281.675016261853\n",
      "Fold: 2. Error: 288.38164186977633\n",
      "Fold: 3. Error: 252.80554916202297\n",
      "Fold: 4. Error: 274.83039775112576\n",
      "Fold: 5. Error: 249.91013817877382\n",
      "LAMBDA AVERAGE: 269.5205486447104\n",
      "Lambda: 0.1\n",
      "Fold: 1. Error: 21.494107244928493\n",
      "Fold: 2. Error: 20.40295905215266\n",
      "Fold: 3. Error: 18.66240844943417\n",
      "Fold: 4. Error: 25.116437906551965\n",
      "Fold: 5. Error: 22.90006062583064\n",
      "LAMBDA AVERAGE: 21.715194655779584\n",
      "Lambda: 0.01\n",
      "Fold: 1. Error: 17.89961806856108\n",
      "Fold: 2. Error: 20.60101543589219\n",
      "Fold: 3. Error: 15.300577288722952\n",
      "Fold: 4. Error: 16.103828700399553\n",
      "Fold: 5. Error: 24.36922875047351\n",
      "LAMBDA AVERAGE: 18.854853648809858\n",
      "Lambda: 0.001\n",
      "Fold: 1. Error: 22.120515998293445\n",
      "Fold: 2. Error: 12.805498902418814\n",
      "Fold: 3. Error: 17.399272285579485\n",
      "Fold: 4. Error: 18.907906539945323\n",
      "Fold: 5. Error: 15.496265314676894\n",
      "LAMBDA AVERAGE: 17.34589180818279\n",
      "Lambda: 0.0001\n",
      "Fold: 1. Error: 18.90513051308685\n",
      "Fold: 2. Error: 17.138574318436756\n",
      "Fold: 3. Error: 24.855574956251054\n",
      "Fold: 4. Error: 18.797927727509116\n",
      "Fold: 5. Error: 18.42840150874861\n",
      "LAMBDA AVERAGE: 19.625121804806476\n",
      "Lambda: 1e-05\n",
      "Fold: 1. Error: 23.92337561024042\n",
      "Fold: 2. Error: 17.214572854892992\n",
      "Fold: 3. Error: 17.61002265462196\n",
      "Fold: 4. Error: 21.151712858887713\n",
      "Fold: 5. Error: 16.60731580859065\n",
      "LAMBDA AVERAGE: 19.301399957446744\n",
      "Lambda: 1e-06\n",
      "Fold: 1. Error: 17.106159456538276\n",
      "Fold: 2. Error: 15.874651443835047\n",
      "Fold: 3. Error: 15.931341544180578\n",
      "Fold: 4. Error: 27.61864882348346\n",
      "Fold: 5. Error: 13.632013387958775\n",
      "LAMBDA AVERAGE: 18.032562931199227\n",
      "\n",
      "\n",
      "**OPTIMAL LAMBDA: 0.001**\n"
     ]
    }
   ],
   "source": [
    "# User must specify lambdas over which to search\n",
    "lambdas = [1, 0.1, 0.01, 0.001, 0.0001, 0.00001, 0.000001]\n",
    "\n",
    "cross_validator = CustomCrossValidator(\n",
    "    X, Y, CustomLinearModel,\n",
    "    loss_function=mean_absolute_percentage_error\n",
    ")\n",
    "cross_validator.cross_validate(lambdas, num_folds=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After identifying the optimal $\\lambda$ for your model/dataset, you will want to fit your final model using this value on the entire training dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([7.75056473, 5.66556818, 8.82981011, 1.08610394, 3.58484621,\n",
       "       5.52614915, 0.21737665, 0.63868509, 3.04998495, 5.39175218])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lambda_star = cross_validator.lambda_star\n",
    "final_model = CustomLinearModel(\n",
    "    loss_function=mean_absolute_percentage_error,\n",
    "    X=X, Y=Y, regularization=lambda_star\n",
    ")\n",
    "final_model.fit()\n",
    "final_model.beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can then generate out-of-sample predictions using this final, fully optimized model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([23.28912996, 30.44766829, 32.25209686, 25.3555125 , 21.03657229,\n",
       "       16.94769912, 20.11264239, 28.44548273, 20.53333071, 28.27208959])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_data = np.random.random((10,10))\n",
    "final_model.predict(test_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
